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ABSTRACT 

We present results from the Spitzer Space Telescope and molecular line observations 
of 9 species toward the dark cloud L43. The Spitzer images and molecular line maps 
suggest it has a starless core and a Class I protostar evolving in the same environment. 
CO depletion is seen in both sources, and DCO"^ lines are stronger toward the starless 
core. With a goal of testing the chemical characteristics from pre- to protostellar stages, 
we adopt an evolutionary chemical model to calculate the molecular abundances and 
compare with our observations. Among the different model parameters we tested, the 
best-fit model suggests a longer total timescale at the pre-protostellar stage, but with 
faster evolution at the later steps with higher densities. 

Subject headings: stars: formation — ISM: individual (L43, RN091) 



1. Introduction 

The first step in our understanding of star formation is to determine the properties of dense 
cores through observations. In general, dust continuum emission is useful to study physical struc- 
tures of the envelope with tools of radiative transfer modeling. Molecular emission lines are comple- 
mentary to trace physical properties through studies of excitation; also the line profiles can uniquely 
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probe kinematics. However, the molecular abundances are not homogeneous for all species through- 
out the core, which easily leads to confusion when interpreting the observed lines. For example, 
the effect of selective depletion is observed: carbon-bearing and sulfur-bearing molecules are eas- 
ily depleted in the inner cold, dense region of starless cores, and nitrogen-bearing molecules can 
remain in the gas phase at higher densities (e.g. Caselli et al. 1999; Tafalla et al. 2002, 2004; 
Lee et al. 2003). Once a central object forms, the chemistry becomes more complicated since the 
temperature distribution can be drastically affected by the additional heating mechanisms such as 
radiation of the central protostellar object and shocks from outflow. The surrounding materials 
are heated up and reach the evaporation temperature of CO, the main gas-phase repository of 
carbon. The abundances of other molecules such as HCO+ and H2CO can be enhanced by their 
desorption off dust grains as well as the higher formation rates due to the increase of gas-phase CO 
abundance. The desorbed CO can also destroy N2H+, resulting in a drop in N2H+ abundance, and 
an anti-correlation of C^^O and N2H"'" emission has been seen toward a few Class protostars (e.g. 
J0rgensen 2004a; Lcc et al. 2007) on 10^ AU scale. Such chemical characteristics obtained from 
molecular spectra can potentially be used to constrain the evolutionary state, but a comprehensive 
picture of coupled dynamical and chemical evolution is needed since the chemistry depends on not 
only elapsed time but also density and temperature. Moreover, our knowledge of critical chemical 
processes, which needs laboratory experiments relevant to interstellar environments, is still limited. 

L43 is a Lynds opacity class 6 (i.e. the highest opacity) dark nebula (Lynds 1962) located near 
the Ophiuchus cloud at a distance of 130 pc (de Geus et al. 1989). A bipolar molecular outflow 
driven by IRAS 16316-1540 (Andre & Montmerle 1994) is seen, as well as a circumstcllar envelope 
elongated perpendicular to the outflow axis (Mathieu et al. 1988; Lee et al. 2002; Lee & Ho 2005). 
The northern outflow lobe is redshifted, and the southern outflow lobe is blueshifted. A reflection 
nebula RN091 is associated with the outflow source (Hodapp 1994). Two sources (hereafter RN091 
and L43E) are detected at submillimeter wavelengths (e.g., Ward-Thompson et al. 1999; Shirley 
et al. 2000; Kirk et al. 2005; Young et al. 2006). It is usually assumed that RN091 is equivalent 
to IRAS 16316-1540, and L43E is a starless core. The mass within 3-o" contours is estimated to 
be 0.3 ± 0.1 Mq for the dense core associated with RN091 and 0.8 ± 0.2 Mq for L43E from the 
Submillimeter Common-User Bolometer Array (SCUBA) 850 //m observations (Young et al. 2006). 
The two sources are roughly 80" apart and provide a test bed for chemical signatures prior to and 
during collapse. 

A few starless cores and cores with central protostars have been previously studied in detail 
with a combination of modeling and observations. Typically, envelope density and temperature 
distributions are determined with dust radiative transfer modeling, which provides the inputs for 
line simulation, and the radial dependence of molecular abundances can be tested. For example, 
empirical abundance models, such as a step function, are used to describe the amount of depletion 
in starless cores (e.g. Lcc et al. 2003). For sources with central heating, a drop function with 
the innermost region undepleted can reasonably reproduce the observed lines (e.g. J0rgensen et 
al. 2004b). Abundance profiles predicted by evolutionary chemical models with varying densities 
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and temperature are also considered; for example, Aikawa et al. (2005) examined the observed 
properties of 4 prcstcllar cores with contracting Bonnor-Ebcrt (BE) spheres, and Evans et al. (2005) 
presented line modeling of a collapsing protostar, B335, with a chemical model incorporating inside- 
out collapse (Shu 1977), described in Lee et al. (2004). In this paper, we also aim at testing such 
dynamical and chemical models, but L43 provides two snapshots, with and without protostellar 
heating, which allow us to compare the chemical behavior at different evolutionary stages based on 
the same initial environmental condition. 

This paper presents new observations of molecular spectra and mid-infrared observations to- 
ward L43, and we follow the method of Lee et al. (2004) in attempts to interpret the data, including 
a dynamical model, continuum radiative transfer, gas energetics, and line radiative transfer. The 
chemical network is updated with a new binding energy of N2 (Obcrg ct al. 2005), and we used 
the N2H+ dissociative recombination rates of Geppert et al. (2004). We have also extended the 
network to include deuterated species since deuterium fractionation is important. This is espe- 
cially true in early phases when it is associated with CO depletion. §2 describes the observations, 
§3 presents observational results and simple analyses, §4 details the modeling sequence, and the 
modeling results are discussed in §5. 



2. Observations 

2.1. Spitzer Space Telescope 

Mid-infrared observations of L43 were performed with the Spitzer Infrared Array Camera 
(IRAC; Fazio et al. 2004) and Multiband Imaging Photometer for Spitzer (MIPS; Rieke et al. 

2004) aboard the Spitzer Space Telescope by the Spitzer legacy Program, Prom Molecular Cores to 
Planet-Forming Disks (c2d; Evans et al. 2003). The IRAC observations (Astronomical Observation 
Request [AOR] ID 0005121280, 0005121792) were obtained in all 4 bands (3.6, 4.5, 5.8, and 8.0 
lim) on 2004 August 16 and 17. Two epochs of observation were taken to allow the identification 
and removal of the asteroids, each with 10" dithering, a 12 second exposure, and a 0.6 second 
of the High-Dynamic Range (HDR) mode image. The MIPS observations (AOR ID 0009409792, 
0009433344) were also taken in 2 epochs at the first 2 bands, 24 /im and 70 /im, on 2004 March 
10. We integrated 1 cycle of 3 seconds for 24 /xm, and 3 cycles of 3 seconds for 70 /xm. 

The data were processed with pipeline version S13.2.0 of Spitzer Science Center (SSC), and 
improved by the c2d team with corrections of artifacts. Detailed discussions of the data processing 
and products by the c2d team is included in the delivery documentation (Evans et al. 2007), as 
weU as in Harvey et al. (2006) on IRAC data and Young et al. (2005) on MIPS data. 

Followup observations that are significantly deeper than those of c2d were performed as part 
of Spitzer program 20386 (PI: Philip C. Myers). The observations with IRAC were made on 
2005 August 22 (AOR key 14612480), and with MIPS at 24 only on 2006 April 8 (AOR key 
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14617600). For the IRAC observations the on-source time of 480 seconds (16 dithers of 30 seconds 
each) is ten times longer that used for c2d. For the MIPS 24 fim observations the on-source time 
is 333 seconds, or about 7 times longer than c2d. 

Followup observations with MIPS at 24, 70 and 160 fim were also performed, as part of Spitzer 
program 30384 (PI: Tyler L. Bourke), on 2007 April 14 (AOR key 18158336). The observations 
utilized scan map mode in all three bands. The data reductions were carried out as described in 
Stutz et al.(2007). The MIPS images are 15'.'5 x 54" in size, with exposure times of ~150 seconds at 
24 fim, ~90 seconds at 70 /im, and ~25 seconds at 160 /im (these are approximate as the sampling 
is non-uniform). The resolutions of the observations are ~ 6", ~ 18", and ~ 40" at 24, 70 and 160 
/im, respectively. 



2.2. Molecular Line Observations 

Line observations were carried out with the CSC10 in the period 1998 July to 2006 June. We 
used a double-side band SIS receiver with an acousto-optical spectrometer (AOS) with a bandwidth 
of 50 MHz over 1024 channels. The pointing uncertainty was about 3.6" on average. The observed 
lines with the reference frequency, the date of observation, the full width at half maximum beam 
size {9h), the main beam efficiency {rimb)^ and the velocity resolution {5y) are summarized in Table 
1. Assuming the source coupling efficiency is unity because the sources are extended, we can 
calculate the radiation (source) temperature (Tr) given by Tr = T^/ij^b, where is the antenna 
temperature corrected for atmospheric attenuation. 

Based on the submillimeter continuum emission mapped with SCUBA (Shirley et al. 2000) at 
450 and 850 fim, we use the (0", 0") position of (16'^34™29.0^ -15°47'3.2") in J2000.0 coordinates, 
and assume that the centroids of RN091 and L43E are at (7", 5") and (89", 6") respectively. The 
observations were performed with position switching at the off position (900", 0"). All spectra were 
then reduced with the analysis routines of the package CLASS. Some of the data were taken with 
different offsets, and in cases with a map, we have synthesized the spectra to the centroid positions 
at a risk of slight degradation of the spatial resolution. 



3. Results 

3.1. Source Properties 

Figure 1 shows a 3-color image of L43 in the IRAC bands 1, 2, and 4. Figure 2 shows the 
850 /xm contours from Young et al. (2006) on IRAC and MIPS images at the 3.6, 8.0, 70, and 
160 /xm bands, as well as a 350 fxm map from Wu et al. (2007) obtained with the Submillimeter 
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High Angular Resolution Camera (SEARCH) at the CSO. L43E is seen as a dark core at 8.0 /um 

and is traced well by the 850 fim dust emission. Photometry including the Spitzer results, the 
submillimeter results, and other observations are listed in Table 2. The table gives the wavelength 
(A), flux density (Si,) with uncertainty (a), and the aperture size. In cases of no detection, an 
upper limit is given. 

RN091 falls in the range of Class I sources, as described in Allen et al. (2004) with the IRAC 
colors [3.6] - [4.5] = 0.93, [5.8] - [8.0] = 0.79. For L43E, there is no point source detected in IRAC 
or MIPS bands. Prom all data of RN091 in Table 2, we also calculate the bolometric luminosity 
(Lfjci), the bolometric temperature (T^o/), and Lgmm/Lboi, where Lgmm includes flux at A > 350jUm. 
The Li)oi is 2.5 it 0.1 L0, and T^oi is 337.6 ± 18.9 K, which fits the criterion for Class I in Chen et al. 
(1995). With L 

smml-^bol — 0.004, RN091 Can also be categorized as Class I as defined by Andre 

et al. (1993). 

3.2. Molecular Lines 

Spectra toward the centroid position of both sources are presented in Figure 3. The line 
parameters derived from the spectra are presented in Table 3, which gives the integrated intensity 
(I), the velocity with respect to the local standard of rest {vlsr), the linewidth (Av), Tr, and the 
number of observed positions. For single-peaked lines, these values are determined from a Gaussian 
fit. For double peaked-lines, Tr is from the stronger peak, v^sR is the velocity of the dip, both 
determined by eye, I is integrated over the full line, and the Av is I divided by Tr. For lines with 
hyperfine structure, vlsr and Av are from hyperfine structure fitting, I is the total area under 
these components, and Tr is for the strongest component. 

The centroid velocities are estimated to be 0.52±0.01 km s'^ for RN091 and 0.77±0.02 km s'^ 
for L43E from the single peaked lines DCO"'" and C^^O, which are likely to be optically thin. Figure 
4 shows the DCO+ map in the center region of L43. The velocities are blueshifted toward the east 
and redshifted toward the west, possibly indicating rotation of the cloud that encompasses both 
sources. The DCO"*" lines are slightly wider near the RN091 peak, but no significant variations in 
linewidths are seen across the mapped region. Figure 5 shows the integrated intensity map of C^''0 
and C^^O; neither of the two CO isotopologues shows line peaks near the center position of RN091, 
which could result from CO depletion, suggesting that the evaporation from central heating is still 
small. Figure 6 shows the 350 fim contours (Wu et al. 2007) overlaid on the extended DCO"^ 
integrated intensity map. Several submillimeter peaks are resolved in the 350 ^m map, and they 
lie on the broad extended structure seen in DCO"*", suggesting that they formed from the same 
structure. CS J = 2 — 1 was mapped by Mathieu et al. (1988) and Lee et al. (2005) for RN091 
tracing both the outflow and the envelope, but the CS J = 7 — 6 line is not detected in our 
observations. 
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3.3. Simple Analysis 

3.3.1. The H2 column density from dust continuum emission 

Assuming the dust grains emit as blackbodies {I^, = Ki,By{TD), B^, is the Planck function), 
and the emission is optically thin at a submillimeter wavelength, we can relate the observed flux 
density {S^,) to the column density of gas (A''(H2)) by 



where ji is the mean molecular weight, mn is the atomic mass, is the dust opacity per gram 
of gas, and Vl is the aperture solid angle. Wc take the 850 /im data in Tabic 2 and measure the 
flux density in an aperture of 34" to match the main beam size of 0^*^© observations. We adopt 
the dust opacities from column (5) of Table 1 in Ossenkopf k, Henning (1994, hereafter 0H5) for 
coagulated dust grains with thin ice mantles, and assume a standard gas-to-dust ratio of 100 to get 
K850 = 0-02 cm^ g~^. Assuming of lOK, we can obtain 

N{R2)rno91,Ss,, = 3-35 X lO^^cm'^; 

iV(H2)L43fi,585o = 4.76 X 10^2 cm-2. 

If Trf is 20K, the calculated N(H2) would decrease by a factor of 3. If we assume a gas-to-dust 
mass ratio of '^124, for example, using the model from Table 4 of Draine (2003; and the references 
therein), the derived iV(i?2) would be 24% higher. 



3. 3. 2. The H2 column density from molecular lines 

Assuming the line is optically thin and all levels in LTE, the integrated intensity for the 
transition J— > J — 1 of a linear molecule x in state and the column density {N(x)) of species x 
is given by 

, SkQeT^ r , 

as in Equation (2) from Lee et al. (2003), where Ej = hBJ{J + 1) in the approximation of rigid- 
rotor, B is the rotational constant, /x is the dipole moment, Q is the partition function, Tex is 
the excitation temperature, and Tr is the measured brightness temperature. We can correct this 
equation for optical depth: 

N(x)thick = N{x)thinZ 

1 — e ^'^ 

We can also convert the column density of a certain molecule into H2 column density by 

Arm ) - 
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where X{x) is the molecular abundance. Under the assumption that all lines are thermalized and 
Tex = WK, we calculate A^(H2) using the observed C^^O J = 2 - 1 line averaged over the beam 
(I in Table 3). We adopt the physical parameters B, fj,, and Q from the JPL catalog, X(CO) of 
2.7 X 10-^ (Lacy et al. 1994), and X(CO)/X(CiSO) ~ 580 (Wilson & Rood 1994). If C^^O is 
optically thin, we obtain 

iV(H2)HiV09i,ciB0 = 1-96 X W^'cm^^; 



cm 



The results are less than one-fifteenth the values calculated from dust emission at 850 //m. 
A lower X(CO), for example, X(CO) ~ 1.8 x 10~^ from Dickman (1978), can result in less CO 
depletion, but the ratio of A^(H2) derived from C^^O to A^(H2) derived from H2 would be still less 
than 0.1 for both sources. Considering that C^^O J = 2 - 1 line might be optically thick, the H2 
column densities derived from the line would be underestimated. To correct the optical depth effect, 
we can fit the C^''0 J = 2 - 1 hyperfine structures to estimate the total optical depth of C^''0 J = 2 
- 1, a nd calculate the optical depth of C^^O J = 2 - 1 with the abundance ratio X(C^^O)/X(C^'^0) 



~ 4 (IWouterloot et al. 20051 ). The hyperfine fitting result of Ten 0(2-1) ~ 0-59 for RN091 has a 
large error due to the low S/N. After the correction for r, A'^(H2) derived from C^'^0 is still lower 
than A^(H2) derived from H2 by a factor of ~ 10, which strongly indicates the depletion of CO 
molecules. Furthermore, the column densities derived from C^^O are similar in both sources, but 
from the 850 /xm data, the column density is higher in L43E than in RN091, which suggests there 
can be more depletion in L43E, with no central heating in the starless core. We can investigate the 
depletion in more detail through the modeling procedure in the next section where the distribution 
of excitation temperature and molecular abundances are taken into consideration. 



3.3.3. Deuterium Enhancement 

At low temperature, the production of H2D^ is favored through the exothermic exchange 
reaction + HD — > H2D"'" + H2. The H2D+ abundance can be further enhanced with CO 
depletion, and its deuteration can propagate to other molecules, such as DCO"*" and N2D+, by ion- 
molecule reactions (Dalgarno & Lepp 1984). Therefore, in cold clouds, the molecular D/H ratio 
can be much higher than the elemental D/H ratio (= 1.5 x 10~^; Oliveira et al. 2003). 

The integrated intensities of DCO^(3-2) toward RN091 and L43E are similar, but the intensity 
of HCO''~(3-2) is much higher in RN091 than in L43E, which must be due to the higher HCO"*" 
column densities and higher gas temperature in RN091. We do not have optically thin lines 
such as H^^CO"*" or HC^^O"*" to help determine the HCO''" column densities (e.g. Appendix Al 
in Caselh et al. [2002]), but we calculate iV(DCO+) and A^(HCO+) for the optically thin case 
(''"z)CO+ (3-2)^1 Thco+{3-2)^^)- For constant abundances and optically thin lines, we obtained 
N{BCO+)l43e/N{RCO+)l43E ~ 4.11 and N{DCO+)rno9i/N(RCO+)rno91 ~ 0.85, assuming 
Tex = 10 K. If Tex is 1 K lower (Tex = 9 K), the column densities can increase by 25% - 30% for 
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each molecule. The optically thin approximation gives a lower limit of 4.8 for the difference of the 
enhancement factor iV(DCO+)/Ar(HCO+) between L43E and RN091. 

In a CO-depleted region, N2D+ can be more enhanced than DCO''" since it suffers from less 

depletion. Assuming the N2H"'" lines arc optically thin, A^(N2H^) toward L43E can be derived to be 
1.23 X W^^cm-^. Also, from §3.3.2, the CO depletion factor (fo) of L43E is ~ 22. From the samples 
of starless cores in Crapsi et al. (2005), objects with higher N2H+ column densities (> 10^'^ cm^^) 
and more CO depletion {/d > 10) tend to have higher A^(N2D+)/A^(N2H+) abundance ratio (> 0.1) 
and can be recognized as the most evolved starless cores. However, the line intensity of N2H+(3- 
2) toward L43E is ~3 times lower than RN091 and N2D+ (3-2) is only marginally detected in 
L43E with a very narrow line. The iV(N2D+)/iV(N2H+) ratio we derived is ~0.03. For RN091, 
the A^(N2D^)/A^(N2H"'") ratio we derived is ~ 0.02, which is at the lower end compared with the 
starless core samples because of heating from the central protostar. 

4. The Modeling Procedure 

In this section we describe the sequence of modeling from the determination of the dust and 
gas temperature distribution to line simulation. Similar techniques were employed in Young et al. 
(2004) and Evans et al. (2005), but in this paper, we apply an updated chemical network. 

For a given density structure (n(r)), we can calculate T^ir) by computing the transport of 
energy through the dust envelope (§4.1), and we simulate the observed spectral energy distribution 
(SED) and the spatial intensity distribution (Ii/(6), where b is the impact parameter). Then Txir) 
can be calculated from a gas energetics code (§4.2). We use an evolutionary model including the 
self-consistent treatment of dynamics, dust continuum radiative transfer, gas energetics, and chem- 
istry (§4.3). Once we obtain the molecular abundances (X(r)), we can calculate the populations 
of relevant energy levels (nj(r)), simulated line profiles, and compare the observations with the 
simulations (§4.4). 

4.1. Dust Models 

We use the radiative transfer package CSdust3 (Egan et al. 1988), which has input parameters 
of luminosity of the central source, dust opacity as a function of wavelength, the strength of the 
interstellar radiation field (ISRF), and the density distribution of the envelope. The program 
ObsSph (Shirley et al. 2002) then takes the output from CSdust3 and simulates the telescope 
observations. For simplicity and because of the limitation of one-dimensional model, we focus on 
the SED fitting at longer wavelengths to trace the envelope. At short wavelengths (A < 60/im), the 
SED is more sensitive to opacity and detailed geometry. For example, RN091 is known to have an 
inclined disk (Weintraiib et al. 1994), and the infrared spectra also show inclination effects, such 
as silicate emission (e.g. Boogert et al. 2008). With an 1-dimensional model envelope, the amount 
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of observed protostellar emission that directly comes through the outflow cavity is not taken into 
account, which can lead to an underestimate of fluxes, particularly at mid-IR wavelengths (e.g. 
Whitney et al. 2003a, 2003b; Robitaille et al. 2006; Crapsi et al. 2008). 

We use the ISRF spectrum based on Cosmic Background Explorer (COBE) data from Black 
(1994) with ultraviolet (UV) wavelengths described in Draine (1978), and modified by dust ex- 
tinction with the properties in Draine & Lee (1984) assuming an attenuation by an external Ay 
outside the core we observe. A scaling factor S[sb,f is also introduced to increase or decrease the 
ISRF at all wavelengths except the cosmic background radiation. For dust opacities in cores, we 
use the model of 0H5 dust since it is successful in reproducing the observed SEDs in star-forming 
cores (e.g. van der Tak et al. 2000; Evans et al. 2001; Shirley et al. 2002). The inner and outer 
boundary of the models are set as R^ = 80 AU and Rout = 8000 AU, which correspond to O'.'G 
and 61'.'5, respectively, and Rout is substantically smaller than the separation of the two sources. 
To compare with models, the observed Iu{b) are calculated from the 850 fj,m maps (Shirley et al. 
2000), and since RN091 and L43E can be partially overlapping, the data in sectors that are possi- 
bly contaminated were excluded. For example, the data in a 90° sector at the east pointing toward 
L43E are cut out in the averaged Iu{b) of RN091. 

Previous observational st udies of pre-protostellar core s (PPCs) suggest a region of roughly 



constant density at small radii (IWard-Thompson et al. 19941 ) , and since they do not have an internal 



source, the only heating source is the ISRF. We attempt to fit the SED and Ii,{b) of L43E with a 
set of BE spheres as discussed in Evans et al. (2001). The central density (ric) of the BE sphere is 
mostly determined by comparing the shape of the normalized luib). The best-fit values are selected 
by minimizing at 850 /xm for luib) as well as Xsed- 

- ^[ ^ ' 

where I°^^{bi) is the normalized intensity at impact parameter 6j, /^°°'(6j) is the modeled intensity, 
(Ti{bi) is the uncertainty of the data, Nb is the number of impact parameters; 

cobs omod 

xIed = .7 f/N^, 

where 5*°^^* is the observed flux, S^°'^ is the modeled flux, as{T^i) is the uncertainty of the data, Ny 
is the number of data points. 

To compute Xsed^ "^^^ i^ excluded because the observing mode of SHARCII 

used for these observations tends to lose flux for extended sources (Wu et al. 2007). The 160 fim 
emission is very sensitive to temperature and is used to constrain the ISRF. We tested models with 
the external Ay in a range between 0.5 to 5.5 and sjsrf between 1 to 10. Figure 7 shows the 
best-fit model with a BE sphere of nc = 3 x lO^cm"^, sjsrf = 5, and the external Ay = 1- Xsed 
for the best-fit model is 7.5 (the blue lines in Figure 7), but different combinations of external Ay 
and sjsRF can provide similar fits, for example, the model with sjsrf = 4 and the external Ay = 
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CS (red lines in Figure 7) and the model with sjsrf = 6 and the external Ay = 1.5 (green lines 
in Figure 7) have Xsed '^bout 15. The external Ay can be further constrained by line modeling 
of C0(2-l) emission since photoelectric (PE) heating is sensitive to the surrounding extinction, as 
detailed in §4.2. 

For the density distribution of RN091, we use the gravitational collapse model presented by 
Shu (1977). In this picture, the model begins from an idealized singular isothermal sphere with an 
envelope with n(r) oc r~^. Collapse begins inside the sphere, and a wave of infall propagates outward 
at the speed of sound as time progresses. The density inside the infall radius {vinf) approximates 
n(r) oc r~^-^ as material falls into the center. For RN091 at A > 60 /xm, Li,oi,\>mnm = 0-8 L©. 
The best-fit model adopting this scenario is shown in Figure 7, with an internal heating source of 
0.4 L0 and rj„/ of 0.006 pc. 

4.2. Gas Energetics 

Deep inside the cloud, the gas temperature is equal to the dust temperature because high 
densities produce rapid collisions between dust and gas. However, at lower densities, since dust 
and gas are decoupled, the gas temperature must be determined by balancing heating and cooling 
rates. We use a gas-dust energetics code by S. Doty including gas-dust energy transfer, cosmic 
ray heating, PE heating, and molecular cooling (details are given in Doty & Neufeld [1997] and 
the appendix in Young ct al. [2004]). CO is the major coolant in clouds, and in principle we 
should include CO depletion in the modeling loop since it could raise Tk in the cold, dense inner 
region. Goldsmith (2001), Doty & Neufeld (1997), and Lee (2004) found that the effect is relatively 
small because of efficient gas-dust coupling at the high densities where CO is significantly depleted; 
therefore, we did not include depletion and desorption in our gas energetics calculation. 

The calculated Txir) is shown in Figure 8. The gas temperature decreases in less dense regions 
where gas-dust collisions are not as efficient, but rises at the outer edge because of PE heating. In 
the gas energetics code, the strength of the UV field relative to the local ISRF on the outside of the 
core is parameterized as Go- We use the equation Go = exp(— 1.8Ay), where the attenuation of the 
ISRF corresponds to a certain external Ay. In dust models, different combinations of external Ay 
and sjsRF can provide similar fits. For example. Figure 7c shows that models with Ay = 0.4 with 
sjsRF = 4, Ay = 1 with sisRF = 5, and Ay = 1.5 with sjsrf = 6 result in similar modeled fluxes 
at 160 /xm. However, the calculated gas temperature proflle is more sensitive to the external Ay 
since the attenuation of the ISRF predominately affects the short wavelengths which are important 
to the PE heating rate. To examine the model temperature at the outer layer and constrain the 
amount of surrounding extinction, we modeled the C0(2-l) line toward L43E and compare with 
the observations (Figure 9), and derive the best-fit Go of 0.165 (corresponds to an external Ay of 
1). Also, the cooling rate depends on the linewidth of CO, and we use a microturbulent width of 
0.24 km s~^ determined by the C^^O line fitting. A more detailed treatment of the line modeling is 
described in §4.4. 
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4.3. Chemical Models 

We use the evolutionary chemical model presented by Lee et al. (2004). A sequence of BE 
spheres with increasing central densities and an inside-out collapse initiated at the point of a singular 
isothermal sphere are combined to describe the dynamics, and the model calculates gas and dust 
temperatures and follows chemical evolution self-consistently at each time step. The collapse model 
is adopted from Young &: Evans (2005) for the evolution of luminosity and density distribution. For 
RN091, we use the physical parameters that best match the dust modeling results in §4.1, and 
of 0.006 pc corresponds to the timestep of 25000 yr. For L43E, since there is no good constraint 
for the timescales and dynamics of the PPC stage, we simply assume 8 steps of different central 
densities and test different timescales for each step to optimize the fitting at ric = 3 x 10^ cm~^. 
In the chemical mo dels, we assume bare Si02 for dust grain surfaces to which molecules freeze 
(|Bergin et al. 19951 ). 

4.4. Line Modeling 

With the best-fit density profile from dust models, the gas temperature profiles from the 
gas-energetics calculation, and molecular abundances {X{r)) from chemical models, the molecular 
excitation was modeled with a Monte Carlo radiative transfer code (mc) which runs with 40 shells 
calculating the level populations. We then simulate the observations with a virtual telescope 
program (vt) (Choi et al. 1995) at the observed offset position. We use the vlsr derived from 
observations for the modeled line velocity in each source. To fit the observed lines, we also use 
simple abundance profiles, such as constant and step function, to compare with the results from 
chemical models. The goodness of a fit is determined by eye. 

We use collision rates by Flower &: Launay (1985) for CO, by Green (1991) for H2CO, by 
Flower (1999) for HCO+, and Daniel et al. (2005) for N2H+. For the deuterated species DCO+, we 
use the de-excitation rate coefficients of HCO"*", and the excitation rate coefficients are obtained by 
detailed balance. In some cases, rates are linearly extrapolated to lower temperatures. For N2II"'" 
and C^''0, the components of the hyperfine structures are partially resolved. Components separated 
by less than the 1/e width of the velocity dispersion were aggregated into a single component. The 
aggregated components are modeled independently with relative intensity-weighted abundances 
and added together to make combined line profiles. 

5. Modeling Results 

The line modeling results are presented with the observations in Figure 3. The abundance 
profiles from chemical models are presented in Figure 10. Some of our line models show self- 
absorption features, which are not seen in the observed line profiles, especially in HCO'*', H2CO, 
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and N2H"'". An explanation for the lack of self-absorption in molecular lines is clumpy cloud 
structure with the existence of macroturbulent motions, which is simulated with the Monte Carlo 
methods for 3-dimensional clouds by Juvela (1997). The line profiles can be confused with heavy 
self-absorption as well as rotation (e.g. Pavlyuchenkov et al. 2008). Since we are not including 
a complete 3-dimensional calculation in this paper, we only fit the line strength at half maximum 
level for lines that are self-absorbed in the models. 

From observations toward star forming cores, one of the features that potentially indicates the 
evolutionary status is an anti-correlation of CO and N2H"'". N2H+ is thought to be a good tracer of 
dense core gas at earlier evolutionary stages until To reaches the temperature for CO evaporation. 
The relative depletion and desorption of N2 and CO are the key processes governing the abundance 
of N2H"'"; however, recent laboratory result shows there is no significant difference of N2 and CO 
binding energies (Oberg et al. 2005; Bisschop et al. 2006). With an updated higher value of N2 
binding energies, the modeled N2H"'" lines are much weaker than the observed lines. The timescale 
for the PPC stage is also a key issue for the freezeout of molecules, so we varied the timescales 
at different central densities. A longer total timescale allows more N2 to form, but a longer time 
at higher densities (n > lO^cm"^) would result in more depletion. N2H+ abundances are directly 
affected by the amount of N2 formation and CO depletion. To fit the observed N2H"'" lines, we 
assume a longer total timescale that allows more N2 formation and CO depletion and shorter time 
at later stages to prevent N2 depletion. Two test models are listed in Table 4. We use the timescales 
in Lee et al. (2004) as the standard model (Model 1), and to optimize the N2H"'" abundances for 
L43E, a longer total timescale of 2 x 10^ yr is assumed (Model 2). In Figure 10, the solid lines are 
abundance profiles from Model 1, and the dashed lines are from Model 2. For both models, we also 
tested the N2H"^ dissociative recombination rates in Geppert et al. (2004 ) to compare with the 



rates in the UMIST database for astrochemistry 1999 (|Le Teuff et al. 20001 ). The rates in Geppert 
et al. (2004) fit the observations better, and are included in our chemical network for the models 
presented in this paper. 

The CO lines are difficult to fit in detail since they are very optically thick and confused 
by the wing components. The line profile toward RN091 has wide outflow wings, and the line 
profile toward L43E also shows a blue wing which possibly comes from another component in the 
foreground along the line of sight. The modeled CO abundances can still be tested through line 
fitting of CO isotopologues. We assume the abundance profiles for the CO isotopologues have the 
same form as CO with an isotope ratio of O/^^O = 580 (Wilson & Rood 1994) and ^^O/^^O = 4 



( Wouterloot et al. 20051 ). In Model 2, the modeled C^^O and C^''0 line intensities are efficiently 
decreased to fit the data since more CO depletion is produced. We also assume an ortho to para 
ratio of 1.5 (Dickens k Irvine 1999) to model the ortho-H2CO (3i,2 - 2i,i) hues. For L43E, the 
H2CO abundance in Model 2 is higher than in Model 1 by a factor of 2, but it is possible that 
self-absorption has greatly weakened the blue peak. 

For the HCO^ lines toward L43E, although the observed line profile might be obscured by 
absorption from a foreground component, the modeled line is still too strong. The cosmic ray 
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ionization rates could affect the results, but neither enhanced nor lowered reaction rates improved 
the fit of both HCO"^ and N2H+ data in both cores simultaneously. We tested empirical abundance 
profiles with lower abundance at large radii. For example, if Gq is higher than wc adopted, the 
higher dissociative recombination rate can lower the abundance in the outer region. However, we did 
not find an improved fit even with extremely low HCO+ abundances in the outer region. To match 
the line intensity according to the red peak, we use a constant abundance X(HCO'*") = 4 x 10~^^, 
which is lower than the chemical model prediction by a factor of ~7.5. For the dcuterated species, 
since they arc favored in a cold region with CO depiction, both DCO"*" and N2D"'" arc enhanced 
relative to HCO+ and N2H+ at high densities. The DCO+ abundances in L43E are increased from 
Model 1 to Model 2 because of the amount of CO depletion, but the modeled DCO+ lines from 
chemical models are still weaker than observed. The DCO+ line emission generated primarily in 
the center region with high densities where HCO+ can be depleted. It is possible that there is 
still more CO depletion than the chemical models predict. A DCO"*" abundance increased by a 
factor of ~8 is needed to fit the observed line intensity for L43E, and for RN091, we use a step 
function with depletion factor of 5 within 0.003 pc. The best-fit empirical abundance profiles are 
shown in Figure 10 with blue dotted lines for L43E and red dotted lines for RN091. To compare 
with the results in §3.3.3, we also calculated the optical depths for the best-fit HCO"^ lines with 
empirical abundances. The optical depth of the simulated HC0''"(3-2) line in RN091 is 23, which 
is much higher than the value of 3 in L43E. This result also supports the analysis in §3.3.3 that the 
A'^(DCO+)/A^(HCO+) factors calculated with optically thin approximation provide a lower limit, 
with iV(DCO+)/A^(HCO+) at least 4.8 times more enhanced in L43E than in RN091. 

To summarize, most of the lines in both sources can be fitted reasonably well with Model 2, 
which has a longer timescale for the PPC stage. However, the change of parameters in our different 
chemical models may only be compensating the effects not yet included. For example, the grain 
properties should also be time-dependent. Grain mantles covered by water-rich ice have strong 
binding energies for molecules to freeze on, but N2 forms slowly and is less affected. As the density 
increases and CO freezes out, the grain mantle surfaces can be contaminated by CO ice. Once grain 
mantles are CO-dominant, the binding energies are lower, which results in less depletion of N2 and 
N2H+ compared with HCO"*", since HCO"*" depletion follows CO depletion. Our modeling results 
of longer timescales at low densities may work equivalently to the effects of higher grain surface 
binding energies (water-dominant) at lower densities. The difference between our modeling results 
and observations may be more specifically explained in a model including time-dependent grain 
properties. Heavy CO (and thus HCO"*") depletion on the initial water-dominant grain surfaces can 
result in more abundant DCO^. With time, grain surfaces are covered with CO (which produces 
CO2 on surface), resulting in less depletion of N2 and more abundant N2H"'". 
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6. Summary 

We have presented Spitzer observations and molecular line observations from the CSO toward 
the two sources RN091 and L43E evolving in the same region. A self-consistent modeling procedure 
is employed to investigate the physical condition in each source, and to compare the chemistry in 
an evolutionary picture. Our results are as follows. 

Prom the Spitzer images, RN091 can be categorized as Class I, and L43E is still starless with 
no embedded source detected. CO depletion is seen in both sources, but the degree of depletion is 
higher in L43E. Also, strong DCO"'" enhancement is seen in L43E. 

Dust modeling results showed that L43E can be fitted by a BE sphere with central density 
of 3 X 10^ cm~^, and RN091 is characterized to have rj„j of 0.006 pc with an inside-out collapse 
model (Shu 1977). For the line modeling, the best-fitting chemical model suggests a longer total 
timescale at earlier stages but fast evolution at higher densities for the PPC stage. The observed 
abundances of N2H+ and CO isotopologues can be matched well in both sources, but in L43E, the 
modeled HCO+ abundance is too high by a factor of 7.5, and the modeled DCO+ abundance is too 
low by a factor of 8. It is possible there is still more CO depletion in the center than the model 
predicted. In this picture, RN091 is possibly the first star forming in this region, and L43E could 
be spending a longer time at the PPC phase. 
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Fig. 1. — Three color image composite of IRACl (blue), IRAC2 (green); IRAC4 (red) toward L43. 
The bright object is RN091, and L43E is the dark region in the east. 



This preprint was prepared with the AAS lATp^X macros v5.2. 




Fig. 2. — SCUBA 850 fim contours on images at (a) IRAC 3.6 /xm from c2d, (b) IRAC 8.0 ^m from 
Spitzer Program 20386, (c) MIPS 70 from Spitzer Program 30384 (d) MIPS 160 nm from Spitzer 
Program 30384, (e) SHARCII 350 fim, and (f) SCUBA 850 fim. 
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Fig. 3. — Molecular line observations toward RN091 and L43E, shown as black solid histograms, 
and the modeled lines from Model 1 (red lines), Model 2 (blue lines), and empirical models for 
HCO"'", DCO"'", and N2H"'" (black dotted lines). Model 2 has a longer timescale for the PPC stage. 
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Fig. 4. — Map of 000"*" (3-2) in the L43 region. The spectra are convolved to the position of RN091 
at (7",5") in 20" steps. The dashed hnes indicate wlsr = 0.77 km s"^ for L43E at (89",6"). The 
velocities toward the east are red-shifted relative to L43E, and the velocities toward the west are 
blue-shifted. 
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Fig. 5. — Integrated intensity map of C-'^^0(2-l) (solid lines) and 0(2-1) (grayscale) convolved 
in 30" steps. The C^^O contour levels start at 0.8 km s~^ and increase by 0.2 km s~^. The C^^O 
contour levels start at 3a and increase by 2a. The star marks the position of RN091, and there 
is no obvious intensity peak in the mapped region in either line. The C^^0(2-l) beam is shown in 
dashed line. 
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Fig. 6. — SEARCH 350/im contours (solid white) overlaid with DCO"'"(3-2) integrated intensity 
map (grayscale). The 350)um contour levels are 2(7, 3a, 5a, and increased by 2a. The DCO"*" 
integrated intensity contours start at 3a and increase by 2a. The DCO+(3-2) beam is shown in 
dashed lines. 




Fig. 7.— Best-fit dust models for botli RN091 and L43E. a) and b): Tlie model fitting of tfie 
normalized /85o(^); c) and d): the SED fitting. The solid blue lines are the model with an external 
Ay = 1 and sjsrf = 5, the dotted red lines are the model with an external Ay = 0.5 and 
sjsRF = 4, and the dashed green lines are the model with an external Ay = 1.5 and sjsRF = 6. 
The different combinations provide similar fits and the modeled lines are on top of each other. The 
short horizontal lines are marks for the upper limits for L43E from Spitzer Observations. 
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Fig. 8. — Density profiles and calculated To and Tk with the same models in Figure 5. a): The 
density profile of BE sphere with ric = 3 x 10^ cm~^; b): the density profile of inside-out collapse 
model with rj^j = 0.006 pc; c) and d): calculated Tjj and Tk- The solid black lines are T£){r) 
from the model with sjsrf = 5 and Ay = 1, the dashed black lines are TD{r) from the model 
with sjsB,F = 4 and Ay = 0.5, the dotted black lines are T£){r) from the model with sisrf = 6 
and Ay = 1-5, the long-dashed blue lines are Txir) from the model with sjsrf = 5 and Ay = 
1, the dot-short-dashed lines are T^ir) from the model with sjsrf = 4 and Ay = 0.5, and the 
dot-long-dashed green lines are Txir) from the model with sjsrf = 4 and Ay = 1.5 
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Fig. 9. — Simple CO line models with the results in figure 7. The blue line is from the model with 
sjsRF = 5 and Ay = 1, the red line is from the model with sjsrf = 4 and Ay = 0.5, and the 
green line is from the model with sjsrf = 4 and Ay = 1.5. X(CO) is assumed 3 x 10~^ in the 
models, but a higher or lower value only affects the size of the self-absorption dip. The modeled 
line velocity is shifted to 0.5 km s~^ to be aligned with the observed line. The blue wing in the 
observed CO line profile might come from another component in the cloud. 



-27- 




-3 -2.5 -2 -3 -2.5 -2 

log r (pc) log r (pc) 



Fig. 10.— Evolution of abundance profiles for a) CO, b) HCO+, c) DCO+, d) H2CO, e) N2, f) 
N2H+, g) N2D"'", and h) at selected timesteps of Model 1 (solid lines) and Model 2 (dashed 
lines). The difference between the two models is that Model 2 has a longer PPC timescale. For t < 
(before collapse), the lines shows BE spheres of central densities 10^ (green), 3 x 10^ (black), and 
10^ (magenta) cm~^. The central densities correspond to different timesteps in each model, and 
the timesteps are listed in Table 4. For t > (after collapse), the timesteps are t = (yellow) 
and t = 2.5 x 10^ (cyan) yr. The blue dotted lines are best-fit empirical models for HCO^, DCO"^, 
and N2H+ in L43E. The red dotted hues are best-fit empirical models for DCO"^ and N2H"^ in 
RN091. 



-28- 



Table 1: CSO Observing Log 



Line 


Frequency 


UT Date 


Beam Width 


riMB 


5y 




(MHz) 




(arcsec) 




(km s ) 


i^CO 2-1 


230537.970 


2005 Mar 


33 


0.82 


0.10 






2006 Jun 


32 


0.77 


0.13 


2-1 


224714.400 


2000 Jun 


33 


0.81 


0.17 




224714.368'" 


2006 Jun 


33 


0.77 


0.13 


C^^O 2-1 


219560.352 


1998 Jul 


28 


0.57 


0.15 






2005 Mar 


34 


0.82 


0.11 






2005 Jun 


34 


0.74 


0.11 






2006 Jun 


34 


0.77 


0.13 


CS 7 - 6 


A cxCiCiCx r\r\r\ 

342882.900 


2000 Jun 


22 


0.77 


0.10 


H2CO 3i2 — 2ii 


225697.787 


2002 Jun 


33 


0.59 


0.15 






2005 Jun 


33 


0.74 


0.13 


HCO+ 3-2 


267557.620 


2000 Jun 


28 


0.81 


0.14 






2006 Jun 


28 


0.77 


0.13 


DCO+ 3-2 


216112.605 


1999 Jul 


35 


0.60 


0.15 






2006 Jun 


35 


0.77 


0.14 


N2H+ 3 -2 


279511.701 


2003 Jul 


27 


0.68 


0.14 


N2D+ 3-2 


231321.635 


2003 Jul 


32 


0.70 


0.13 






2005 Jun 


32 


0.74 


0.11 



"Reference frequency for the hyperfine shifts in Ladd et al. (1998) 
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Table 2: Observed Flux Densities 



Source 


A 




Aperture 


Reference 




(^m) 


(mJy) 


(arcsec) 




RN091 


1.6 


250(10) 


12 


1 




1.65 


143(6.86) 


2.5*^ 


2 




2.17 


274(8.34) 


2.5*^ 


2 




2.2 


480(20) 


12 


1 




3.4 


510(20) 


12 


1 




3.6 


737(51.1) 


1.7'^ 


3 




4.5 


1110(66.9) 


1.7*^ 


3 




5.8 


1440(72.0) 


1.9'^ 


3 




8.0 


1660(82.8) 


2.0*^ 


3 




24 


5000(200)^= 


70 


4 




70 


20300(4060)'= 


70 


4 




160 


67400(13480)^= 


96 


4 




190 


38000(15200) 


60 


5 




350 


9900(1600) 


40 


6 




450 


11800(2000) 


40 


7 




850 


1600(90) 


40 


7 


L43E 


3.6 


< 6.99^ 


1.7^ 


3 




4.5 


< 10.17<i 


1.7^ 


3 




5.8 


< 57.3^^ 


1.9'' 


3 




8.0 


< 61.5<i 


2.0*^ 


3 




24 


< 378^ 


5.7^ 


3 




70 


400(80)^= 


70 


4 




160 


37000(7400)^= 


128 


4 




350 


6500(1300) 


40 


6 




450 


11400(1900) 


40 


7 




850 


1800(100) 


40 


7 



"FWHM of 2 Micron All Sky Survey point-spread profile 
''FWHM of Spitzer point-spread profile 

■^Uncertainty given in c2d delivery documentation (Table 21) 

'*3-cr value estimated with the Sensitivity Performance Estimation Tool (SENS-PET) from Spitzer Science Center 

REFERENCES.— (1) Myers ct al. 1987; (2) 2 Micron All Sky Survey; (3) results from Spitzer c2d; (4) results from 
Spitzer Program 30384; (5) Ladd et al. 1991; (6) Wu et al. 2007; (7) Shirley et al. 2000 
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Table 3: CSO Obsorvalioiial Rcsullts 



kjUUl Ct; 


Line 




vlsr 


i-\ V 


-L R 


UliocU 








I XV Kill } 


1 Kill t> 1 


{ l.^'Tn t? ^ \ 
1 Kill J 1 




1 T' C *3 I 








99 i^n^'n 1 Q^i 




9 P.Q(C\ r\A\ 
z.oy i^u.uo j 


7 S7('n 91 


(^0, 6} 


1 






u.oy I u.uo 1 


U.O 1 I U.UZ I 


U.^ol U.UU J 


u.^y I u.uo J 




1 fiQ 
luy 




pisn 9 - 1 


1 71 Dfi"! 
±. 1 11 U.UU 1 


n ^1 01 


u.oy I u.uo J 


9 7A(r\ 1 P,^ 


by llLllCblZfcJQ lU (l, OJ 


1 91 




pq 7 _ fi 








\U.^U ^1 .lll.b U.ll) 


10, oj 


1 




il2^VJ 012 ^11 




U.OO^U.UO ) 




1 91 (r\ 97"! 




1 




1 1 V V / o — z 




U.4o(^U.Ui j 








1 




nr!o+ 3-9 








9 q9fO 1 6^1 


Q'^T'Ti't' npQi^firi (7 ^1 

O y llLllCOlZyCU tVJ v ' ) / 


342 




N2H+ 3 -2^= 


3.78(0.07) 


0.40(0.01) 


0.64(0.03) 


3.29(0.07) 


(7, 5) 


11 




N2D+ 3-2 


0.21(0.03) 


0.18(0.03) 


0.37(0.06) 


0.53(0.01) 


(7, 5) 


5 


L43E 


12C0 2 1 


10.10(0.13) 


0.46(0.01) 


1.38(0.02) 


6.88(0.15) 


(116, -13) 


5 




Ci*^0 2 - 1 


1.40(0.09) 


0.66(0.02) 


0.55(0.05) 


2.38(0.18) 


(96, -13) 


6 




H2CO 3i2 - 2ii 


0.47(0.04) 


0.59(0.03) 


0.58(0.07) 


0.76(0.01) 


(89, 6) 


1 




HCO+ 3-2 


1.11(0.07) 


0.94(0.01) 


0.47(0.04) 


2.24(0.10) 


(89, 6) 


5 




DCO+ 3-2 


1.65(0.05) 


0.88(0.01) 


0.50(0.02) 


3.13(0.10) 


(90, 10) 


1 




N2H+ 3-2"= 


1.19(0.10) 


0.73(0.01) 


0.23(0.03) 


1.36(0.10) 


(86, 3) 


25 




N2D+ 3-2 


0.11(0.03) 


0.48(0.03) 


0.26(0.07) 


0.39(0.11) 


(96, -13) 


5 



"(0,0) is the reference position of (16''34'"29.0" , -15°47'3.2") (J2000.0), and the centroids of RN091 and L43E are at 
(7, 5) and (89, 6), respectively. 
''Number of observed positions. 

'^Double-peaked line: Tr refers to the stronger pealc, ylsr refers to the velocity of the dip, I is of all area, and Av is 
I divided by I.QQTr for the gaussian equivalent linewidth. 

'^Hyperfine structure: vlsr and Av are from hyperfine structure fitting, I is of all area, and Tr refers the main peak. 



Tal)le 4: Chemical Models 





iVloaoi i 


iVlouel z 


Central Density 


Tinicstep 


Timestep 


(cm-3) 


(yr) 


(yr) 


10^ 


5 X 10^ 


3 X 10^ 


3 X 10'^ 


1 X 10"* 


5 X 103 


10*^ 


2 X 10"* 


1.2 X 10* 


3 X 10^ 


3 X 10* 


3 X 10* 


10^ 


6 X 10* 


10^ 


3 X 10" 


1.25 X 10^ 


3.5 X 10^ 


10* 


2.5 X 10^ 


1.5 X 10^ 


103a 


2.5 X 10^ 


1.5 X lO'^ 



"The initial constant density. 



